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1.  Introduction 


As  numerical  models  of  the  atmosphere  become  more 
finely  resolved,  restraints  on  the  size  of  time  steps  become 
increasingly  difficult  to  deal  with.  These  restraints  in 
explicit  time  differencing  are  due  to  the  well-known 
Courant-Friedrich-Lewy  (CFL)  conditions  which  limit  the  time 
step  to  the  ratio  of  the  spatial  increment  to  the  phase 
velocity  of  the  predicted  variable.  Thus,  as  the  spatial 
increment  is  decreased,  the  time-step  must  also  be  reduced. 
Implicit  time  schemes  avoid  these  restrictions.  It  has  been 
shown  that  for  linear  systems,  time  schemes  where  the  time 
dependent  variable  in  the  numerical  approximation  to  the 
analytic  eguation  is  given  at  the  predicted  time-step 
(thereby  requiring  an  "implicit"  solution)  are  not  subject 
to  CFL  limitations.  In  fact,  an  Euler-backward  scheme  where 
all  references  to  the  variable  are  given  at  the  predicted 
time  step  can  be  shown  to  be  absolutely  stable  no  matter  how 
large  the  time  increment  relative  to  the  space  increment. 
This  characteristic  makes  implicit  schemes  appealing  for 
application  to  numerical  weather  prediction.  But  there  are 
other  factors  involved  in  numerical  weather  prediction.  For 
one  thing,  the  equations  are  non-linear  and  no  one  can 
guarantee  stability,  even  with  implicit  schemes,  for  non¬ 
linear  formulations.  Non-linearity  also  poses  practical 
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problems  in  solving  for  the  predicted  time  step.  Whereas 
with  linear  equations,  the  terms  involving  the  predicted 
variable  can  easily  be  separated  from  terms  involving  known 
quantities,  with  non-linear  configurations,  the  separation 
and  hence  the  solution  methodology  may  not  be  so 
straightforward.  Some  of  these  impediments  can  be  overcome 
by  introducing  semi-implicit  schemes. 

Semi-implicit  schemes  derive  their  label  from  the  fact 
that  only  portions  of  the  equation  (generally,  the  linear 
portions)  are  cast  in  implicit  formulations.  The  rest  of 
the  equation  is  left  as  an  explicit  expression  of  known 
quantities.  In  this  form  the  scheme  is  not  absolutely 
stable,  and  the  size  of  the  time-step  does  play  a  role. 
Kurihara  (1965),  in  fact,  demonstrated  that  semi-implicit 
schemes  are  weakly  unstable  and  may  not  be  able  to  resolve 
short  waves  efficiently,  if  the  time  step  is  excessive.  The 
problem  is  magnified  when  physical  processes  are  included. 

In  choosing  a  semi-implicit  scheme  for  the  relocatable 
limited-area  model  (RLAM)  developed  by  ST  Systems,  Inc. ,  as 
described  by  Tung  et  al.,  (1988),  in  support  of  the  global 
modelling  effort  at  the  Air  Force  Geophysics  Laboratory 
( AFGL) ,  several  different  models  were  examined.  These 
included  AFGL's  global  spectral  model  (gsm)  itself,  the 
baroclinic  model  of  the  Met  Service  of  Canada  (Robert,  et 
al.,  1972),  and  the  regional  model  of  the  Australian 
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Numerical  Meteorological  Research  Centre,  now  Bureau  of 
Meteorology  Research  Centre  (BMRC)  as  detailed  by  McGregor, 
et  al.,  (1978).  The  methodology  in  all  three  models  is  the 
same.  A  single  equation  in  one  variable  is  solved.  The 
solution  is  then  used  to  solve  for  the  other  variables 
expressed  in  terms  of  the  first.  In  the  case  of  the  gsm,  an 
updated  value  of  divergence  is  derived  first;  afterwards 
vorticity,  surface  pressure,  and  temperature  are  forecast 
based  on  the  new  value  of  divergence.  The  Canadian  model 
first  solves  for  a  variable  which  combines  geopotential  and 
surface  pressure  and  derives  all  others  from  it,  while  the 
Australians  first  solve  for  temperature.  In  all  cases,  the 
single  equation  in  one  variable  is  a  Helmholtz  equation, 
where  the  right  hand  side  contains  all  the  non-linear  terms 
evaluated  from  known  quantities  at  the  present  or  past  time 
steps.  The  gsm  and  the  Canadian  model  approach  were 
rejected  in  favor  of  the  Australian  model  because  of  greater 
compatibility  with  RLAM.  The  gsm,  being  a  spectral  model, 
can  carry  vorticity  and  divergence  as  principal  dependent 
variables  that  are  solved  by  the  model's  algorithm  and  then 
transformed  diagnostically  into  wind  velocity  components. 

The  RLAM,  being  a  primitive  equations,  finite-difference 
model,  as  shall  be  seen  in  the  next  section,  carries  the 
velocity  components  themselves  as  principal  variables.  The 
Canadian  model  does  employ  finite  differencing,  but  its 
vertical  staggering  of  variables  is  dissimilar  to  the  gsm 
and  RLAM.  It  may  have  been  possible  to  reconfigure  the 


Canadian  model  for  use  by  RLAM,  but  since  BMRC's  model  was 
completely  compatible  with  RLAM,  the  path  of  least 
resistance  was  chosen.  The  construction  of  the  BMRC  model 
is  theoretically  consistent  with  RLAM ' s  and  and  as  such,  it 
was  expected  that  adapting  to  RLAM  would  pose  no  problem. 
Unfortunately,  as  one  discovers  quite  often  in  numerical 
modeling,  theoretical  consistency  does  not  always  lead  to 
applicability.  In  the  ensuing  sections,  difficulties  in 
adapting  the  BMRC  scheme  to  RLAM  will  be  presented. 

2 .  Review  of  RLAM 

RLAM  has  been  documented  in  various  stages  of  its 
development,  notably  in  Tung  et  al.,  (1988)  and  Halberstam 
(1988).  Briefly,  RLAM  was  designed  to  have  multiple  modules 
so  that  forecasts  can  be  made  using  any  combination  from  a 
number  of  differencing  schemes  both  spatial  and  temporal. 

The  spatial  choices  consist  of  second  or  fourth  order 
differencing.  Fourth  order  compact  differencing  had 
originally  been  included  but  has  been  dropped  because  it  did 
not  offer  much  in  the  way  of  advantage  over  the  regular 
fourth  order  scheme  and  required  specifications  along  the 
boundaries  which  were  just  as  troublesome  to  deal  with  as 
the  regular  fourth  order.  Along  with  the  horizontal 
differencing,  one  can  specify  the  frequency  and  order  of 
smoothing  along  the  lateral  boundaries  and/or  over  the 
entire  field.  Second  or  fourth  order  diffusion  can  also  be 
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specified  with  a  variable  (in  space)  diffusion  coefficient. 
Physical  parameterizations  are  optional.  The  model  allows 
invoking  one  or  all  of  bulk  boundary  layer  specifications, 
large-scale  precipitation,  and  a  Kuo-type  convection. 

The  lateral  boundary  conditions  are  also  subject  to  the 
selective  process.  At  present,  one  can  select  among  a 
Davies-type  (1976),  Perkey-Kreitsberg  (1976),  or  Orlanski- 
type  (1976)  boundary  specification.  The  selection  of 
lateral  boundary  specifications  is  crucial  for  limited-area 
models.  Great  care  must  be  observed  in  choosing  boundary 
conditions  which  are  compatible  with  the  spatial  and 
temporal  differencing  schemes.  In  any  case,  because  gsm 
data  form  some  part  of  the  boundary  specification,  it  is  not 
conceivable  that  the  limited-area  model  can  depart  too  far 
from  the  gsm  without  causing  major  disruptions.  Thus,  one 
is  faced  with  the  disparate  goals  of  expecting  the  limited- 
area  model  to  improve  over  the  global  model  in  the  interior 
on  the  one  hand  while  also  expecting  the  global  model  to 
provide  useful  boundary  conditions  for  the  limited-area 
model.  Because  in  practice  the  model  forecasts  tend  to 
drift  apart,  modelers  have  had  to  devise  boundary 
specifications  of  the  form  mentioned  here,  where  either  a 
slow  blending  of  large-scale  and  interior  boundary 
conditions  occurs  or  internal  waves  are  allowed  to  exit  the 
domain  without  much  disruption  while  exterior  waves  are 
allowed  to  enter  if  conditions  of  this  sort  prevail.  It 
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also  requires  that  the  forecast  be  limited  in  time  so  that 
the  boundaries  do  not  swamp  the  interior  completely. 

For  time  schemes,  RLAM  originally  included  either  a 
standard  explicit  leap  frog  scheme  or  a  Brown-Caropana 
pressure-gradient  averaging  scheme.  The  latter  is  still 
explicit  but  the  time-averaging  affords  more  liberal  time 
steps.  Because  these  time  schemes  are  explicit,  the  time 
steps  have  to  be  shortened  as  the  spatial  resolution  is 
increased.  Thus,  for  high-resolution  integrations,  a  sern^.- 
implicit  or  equivalent  scheme  was  required  in  order  to  keep 
the  model  tenable  and  economically  feasible.  As  already 
mentioned,  the  BMRC  scheme  was  selected  because  c  its 
compatibility  with  RLAM.  The  following  section  summarizes 
the  theoretical  basis  for  the  scheme. 


3.  The  BMRC  semi-implicit  scheme 

Starting  with  a  set  of  primitive  equations,  the  BMRC 
model  consolidates  the  whole  set  into  a  single  equation  in 
one  variable.  The  original  equations  are 
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T  -  temperature 

u  -  Eastward  wind  velocity  component 
v  -  Northward  wind  velocity  component 
q  -  specific  humidity 
p.  -  surface  pressure 
0  -  geopotential 
t  -  time 

x,  y  -  horizontal  map  coordinates 

p 

a  -  vertical  coordinate  =  p-  (p  is  pressure) 
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f  -  coriolis  parameter 

mx ,  -  horizontal  map  factors 

P  K  R 

7T  -  Exner  function  (g-)  ;k=<£- 

Fo  P 

R  -  Ideal  gas  constant 

cp  -  heat  capacity  of  air  at  constant  pressure 

a  -  vertical  velocity  »  ^ 

Fx ,  F  -  horizontal  friction 
H  -  diabatic  heat  sources  or  sinks 
S  -  moisture  sources  or  sinks 
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Using  the  following  notation  for  time  averaging,  — - —  =  a  , 

allows  us  to  write  the  discrete  approximation  to  the  time 
derivative  as  (A*  -  At-1)/At,  where  A1*1  is  the 

variable  A  at  time  t+At  and  A*'1  is  A  at  time  t-At.  This 
centered  differencing  along  with  the  accompanying  right-hand 
side  containing  terms  in  A1,  requires  specification  at 
three  time  levels.  We  also  separate  the  temperature,  T, 
into  a  climatological  component,  T0,  which  depends  only  on 
the  vertical  coordinate,  a,  and  a  spatially  and  temporally 
varying  component  T' .  Given  these  specifications  we  can 
approximate  (1)  as 
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or,  m  finite  difference  specification.  At  W.  =  In  p_  -  In  p;  • 
The  moisture  equation  is  solved  explicitly  using  the  updated 
values  of  u  and  v  for  advection  and  is  therefore  not 
involved  in  these  manipulations  and  hence  not  listed  in  (2). 
a,  b,  and  c  account  for  the  non-linear  portions  of  the 
momentum  and  thermodynamic  equations  and  are  given  as 
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where  variables  without  superscripts  are  to  be  evaluated  at 
time  t.  We  must  still  specify  a  vertical  differencing 
scheme;  horizontal  differencing  can  be  second  or  fourth 
order  as  desired.  If  we  stagger  the  vertical  so  that  w  is 
specified  at  the  a  levels  or  interfaces  (to  be  designated  a) 
as  opposed  to  all  other  variables  which  are  carried  at  the  a 
layers,  the  vertical  derivatives,  except  possibly  those  of 
w,  should  have  obvious  finite  depictions. 

If  the  divergence  (i.e.,  the  second  term  on  the  left- 
hand  side)  in  (2)d  is  replaced  by  the  differentiated  linear 
pressure  gradient  terms  on  the  left-hand  sides  and  the  non¬ 
linear  terms  on  the  right-hand  sides  of  (2)a  and  (2)b,  we 
arrive  at  the  following  relationship; 
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The  last  term  on  the  left-hand  side  comes  from  a  combination 
of  the  second  derivatives  of  In  p_  at  times  t  and  t-At  and 
substituting  the  aforementioned  finite-difference 

— r  — t  _T 

approximation  to  W.,  i.e.,  AtW.  =  In  p.  -  In  p*  . 

The  V2  operator  is  a  generalized  operator  and  is  given  by 
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where  A  is  any  variable.  Note  that  when  the  mapping  is 
conformal,  i.e.,  m*  =  niy,  the  operator  reduces  to  its  more 
familiar  two-term  Laplacian. 

After  choosing  a  horizontal  differencing  scheme, 
vertical  finite  differencing  will  determine  the  discretized 
form  of  (4).  Remembering  that  W  is  staggered  in  the 
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vertical  with  respect  to  the  other  variables,  m  layer 
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k  can  be  discretized  by  (Wk+1  -  Wk)(Aak)'1.  (4)  is  then 

represented  by  this  matrix  form: 

C(WJ  -  -  E  ?2(WJ  =  -m*Vdk>  (5) 

where  the  block  letters  are  KxK  matrices  K  being  the  total 
number  of  layers,  and  (k)  represents  a  Kxl  vector.  Here 


< 


and  dk  represents  the  right-hand  side  of  (4)  divided  by 
-n^niy.  Note  that  although  there  are  K+l  levels,  (Wk)  is 
only  a  vector  with  K  components  because  Wk+1  =  0  and  need 
not  be  included. 

( 2 } c  can  now  render  a  relationship  between  (Tk)  and 
(Wk)  ,  while  the  hydrostatic  equation  relates  (Tk)  to  (*k)  . 
Vertical  differencing  of  (2)c  can  be  accomplished  in  a 
number  of  ways,  however,  and  the  choice  bears  significant 
consequence.  Because  T0  is  specified  arbitrarily,  one  can 
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define  T0  and  analytically  at  both  layers  and 
levels.  Thus,  while  (2)c  must  be  specified  at  a  layer  and 
therefore  requires  some  kind  of  averaging  of  the  level- 
defined  W,  the  coefficients  of  W  can  either  be  defined 
directly  at  the  layers  or  defined  at  the  levels  and  averaged 
along  with  W.  Instinctively,  it  was  felt  that  the 
differencing  on  the  left-hand  side  of  the  equation  should  be 
consistent  with  the  differencing  employed  for  the  right-hand 
side  (rhs) ,  given  in  (3)c.  But  here  again  a  number  of 
choices  are  available  even  for  the  rhs.  Since  W  is  defined 
at  the  levels  while  T'  is  defined  at  the  layers,  we  can 
average  W  over  surrounding  levels  and  use  centered 
.  3T ' 

differencing  for  .  This,  in  effect,  is  the  method 

employed  by  BMRC  (not  mentioned  in  the  1978  article) .  But 
this  can  lead  to  problems  at  the  first  or  uppermost  layer, 
where  T'  is  not  defined  below  the  first  or  above  the  top 
layer.  BMRC,  in  fact,  extrapolates  downward  to  estimate  a 
surface  temperature  and  assumes  a  constant  lapse  rate  at  the 
top.  One  can  avoid  this  questionable  procedure  by  simply 
approximating  the  terms, 


3T ' 


da 


w 
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for  any  layer  k  as 


|(x  T^HW^  +  Wk)  -  §[(Tk,1-Tk)Wk+,<ok„  -  ok)  ■’  + 

(Tk  —  Tk .  1  )  Wk  ( C7k  <7k .  I  )  ]  +  2  £  ^  1  ”  )  ^k«1  (°k*1  ~  °k  )  + 

(Tk  -  T'k.,)ak(ak  -  <?(,.,)  1  ]W. 

Now,  since  Wk4l  and  aK4l  are  o,  there  is  no  need  for 
defining  a  TK+1.  On  the  other  hand,  a,  =  1  and  W,  =  W.  so 
that  the  second  of  the  bracketed  terms  will  cancel  when  k=l 
and  there  is  no  need  to  define  T0.  Representation  of  the 
first  term,  given  analytically  as  kT'o’w,  is  also  possible 
in  more  than  one  way  because  a  can  be  defined  either  at  the 
level  or  layers.  Thus  the  first  term  may  also  be  written 

i  t  w  W  \ 

i(K  Tk)  +  -r*  As  we  shall  see,  the  choice  of  this 

*  lak*i  <V 

representation  can  have  a  profound  effect  on  the  results  of 
the  model. 

With  any  of  the  discretization  specified  above,  (2)c 
can  be  written  in  generalized  matrix  form  as 

W  -  A(WJ  =  (ck)  (6) 

where  the  matrix  A  is  defined  as  a  KxK  matrix  with 
elements 
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sf  0  0  0  .  0 

s]  sf  0  0  .  .  0 

s3  SJ  0  .  •  •  o 


where,  for  example,  if  we  take  the  non-centered  form  over 
the  layers,  then 

Gk  =  ^k»1  (Tok  +  1  ~  (CTk*1-<7k)  +  ^k^ok  ~  Tok.1)(°k  “  °k- 1  ) 

i  kT„v 

sl  =  -of  -  (Tok  -  T0k  ,)(<7k  -  ok.,) 

S k  Ai—  -(Tok*i  -  Tok)(ok0  -  ok)  1 

The  hydrostatic  equation  in  matrix  form  can  be  represented 
as 


A 


At 

2 


Gl  +  Si 

G  , 


GK 


(*k  -  *.)  =  B(Tk)  (7) 

where  the  matrix  3  is  dependent  on  the  particular 
discretized  relationship  one  chooses  between  layer  heights 

— r  _t 

and  temperature.  One  can  now  substitute  for  (*k)  and  (Wk) 
in  (5)  using  (6)  and  (7) ,  and  arrive  at  one  equation  for 
(Tk)  ,  i.e., 
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(B  +  EA’V(tI)  -  CA"1^)  =  mKmy(dlt)  +  EA  V(ck)  -  CA  ^c,)  (8) 


By  multiplying  all  terms  by  (t  +  EA'  )  1  and  combining  all 
terms  on  the  right-hand-side,  one  remains  with  a  classical 
Helmholtz  equation  in  three-dimensions, 

?2(TTk)  -  GW  =  <Rk)  (8a) 

Solution  of  (8a)  first  requires  reduction  to  a  two- 
dimensional  problem  by  projecting  the  vectors  onto  each  of 

the  eigenvectors  of  the  matrix  Q.  This  is 

—  T 

accomplished  by  substituting  for  (Tk)  and  (Rk) 
by  (rk)  and  (Tk)  ,  where  the  latter  are  the  former 

times  the  inverse  matrix  of  eigenvectors  of  i.e., 

(Tk)  =  ]}{\)> 

where  \  is  the  matrix  whose  columns  are  the  eigenvectors  of 
(}.  (8a)  then  becomes,  after  multiplying  by  V  , 

V2(rk)  -  A(*k)  =  (rk), 

where  is  the  diagonal  matrix  of  eigenvalues 
corresponding  to  the  eigenvectors.  This  equation  neatly 
divides  into  K  two-dimensional  Helmholtz  equations  which  can 
normally  be  solved  by  both  direct  and  indirect  methods. 
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However,  because  the  map  factors  appearing  in  the  Laplacian 
operator  are  spatially  dependent,  only  indirect  iterative 
methods  of  solution  can  be  invoked. 

.  .  — T 

After  t  is  determined,  T  is  found  by  multiplying 
through  by  the  eigenvectors.  Once  T*  is  known,  one  can 
theoretically  invert  (6)  to  derive  (Wk)  and  hence  In  p*, 

derive  (*k)  from  (7)  and,  finally,  u  and  v  from  (2) a  and 
(2)b.  This  is  the  procedure  suggested  by  BMRC,  but  there 
seem  to  be  some  latent  problems  associated  with  this 
algorithm,  as  will  be  described  in  the  next  section. 

4.  Deriving  pressure  from  temperature:  some  pitfalls 

— T  _T 

To  derive  W  from  T  one  inverts  (6)  to  obtain 

(wk)  =  aN*;  -  ck)  o) 

The  updated  surface  pressure  can  then  be  derived  by  invoking 
the  finite-difference  approximation  mentioned  earlier,  i.e., 

AtW*  =  In  pT  -  In  p*'1. 

•  • 

At  this  point,  we  may  consider  what  variables  have  yet 
to  be  solved  in  the  system  of  equations  (2)  before  stepping 
forward  to  the  next  time  period.  Methods  for  deriving 
u**1,  and  v1*1  have  been  outlined  in  the  previous  section. 

But  the  vertical  advection  terms  in  (3) A  and  (3)B  are  given 
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in  terms  of  a,  which  has  not  been  explicitly  specified.  In 
truth,  a  was  preempted  by  the  creation  of  W  to  serve  as  a 
substitute  variable.  This  exchange  should  require  one  to 

replace  a  in  the  non-linear  terms  with  W  similar  to  (3)C. 

If  a  is  kept,  one  may  regard  the  equation  defining  W  in 
terms  of  a  as  a  diagnostic  equation  for  b.  If  so,  (9)  is 
ipso  facto  a  solution  for  6,  as  well.  Yet  BMRC  rederives  a 
from  the  continuity  equation  using  updated  values  of  u  and 
v.  This  unnecessary  computation  will  not  always  yield  the 

same  a  values  as  would  be  derived  from  W  and  is  not 
consistent  with  the  rest  of  the  model.  However,  when  we 

attempted  to  derive  a  directly  from  W,  we  found  the  values 
to  be  unreasonably  large,  as  was  the  value  of  W.  which  is 
the  surface  pressure  tendency. 

What  caused  these  problems  seems  to  stem  from  the 
inversion  of  the  matrix  f[  in  (9)  to  obtain  (Wk)  .  For  the 
sake  of  demonstration,  we  may  choose  a  constant  T0  that  does 
not  vary  with  height  without  any  loss  of  generality.  In 

fact,  the  nature  of  /\  and  its  inverse  depends  little  on  T0> 
as  one  would  expect.  The  values  of  Gk  in  the  definition  of 
A  then  all  become  zero,  and  fa  becomes  a  bidiagonal  matrix 

with  elements  sk  and  sj|  equalling  cr' 1  or  ak\  depending 
on  the  choice  of  discretization.  One  would  also  hope  that 
it  should  matter  little  whether  we  discretize  the  term 

k  T0  o'1  as  icT^1 ,  as  T,,^1  +  ok«i)#  or  as  anY 

other  weighted  combination  of  bk  and  ak.].  But  the 
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If  we 


inverse  of  ^  is  profoundly  affected  by  this  choice. 


select  the  second  option,  i.e.,  k  T0  Wo 
then 


-  1  'N'  1 

=  To 


w 


k  + 1 


'k«1 


2 

AtxT0 


a 


-a 


2 


±b2 


0  0  0  0 


iff 


10 


with  the  last  column's  sign  dependent  on  whether  K  is  odd  or 

—  T 

even.  W.  depends  only  on  the  first  row,  which,  because  a, 

=  1,  consists  merely  of  alternate  entries  of  ±1.  Thus,  the 
pressure  tendency  will  depend  exclusively  on  the  magnitude 

—  T 

and  pairing  of  the  differences  between  T  and  c,  the  non¬ 
linear  tendency  of  the  temperature  including  physical 
processes.  For  pressure  tendencies  to  remain  reasonable, 

one  must  require  that  2AtW,  be  of  the  order  10  .  Since 

-7  .  -?  — T 

(kT0)  is  of  the  order  of  10  ,  T  -  c  should  be  close  to 

or  less  than  10  C  in  order  for  the  pressure  tendency  to 
remain  reasonable.  Note  that  there  is  no  apparent 

—  T 

dependency  of  surface  pressure  tendency  on  At  as  2AtW. 

l  i  ■ 1 

cancels  the  ^  in  the  coefficient  of  f\  .  However,  At  is  a 
factor  in  c,  as  given  by  (3)C,  underlining  the  fact  that  the 
time  differencing  is  not  fully  implicit,  so  that  the  time 
step  still  affects  the  non-linear  tendency  but  not  the 
linear  tendency.  Because  c  contains  physical  processes,  it 
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is  unreasonable  to  expect  that  T  -c  be  constrained  to  be 

less  than  102°C  at  all  times  and  it  is  certainly 
unreasonable  to  expect  that  there  be  any  specific  coupling 

— T 

between  layers  so  that  the  effects  of  T  -c  cancel  each 
other.  This  is  especially  true  near  the  surface  where 
condensation  and  boundary  layer  processes  affect  the  first 
layer  temperature  more  than  other  layers. 

If  the  first  finite  differencing  is  invoked,  i.e., 

T0  +  Wk),  then  A"  becomes 

|  a,  -a2  a ,  -a, 

__2 _ '  0  o2  o3  a, 

KToAt 

0  0  0  0 

Here,  as  the  a  values  decrease  with  increasing  k,  the 

—  T 

stability  can  be  maintained  even  with  higher  values  of  T  -c 
at  upper  layers.  This  would  point  to  a  vertical 
discretization  with  higher  resolution  near  the  bottom  and 
less  resolution  near  the  top  of  the  model.  Thus,  by  what 
seems  like  a  simple  assumption  with  respect  to  vertical 
finite  differencing,  appreciable  consequences  with  regard  to 
stability  can  occur.  This  would  point  to  an  extreme 

sensitivity  on  the  part  of  the  matrix  This  sensitivity 
suggests  an  ill-conditioned  situation  where  small 
differences  in  data  are  magnified  to  levels  far  beyond  what 
one  would  expect. 
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As  a  simple  example  consider  this  12-layer  profile  of 

— T 

T  -c  taken  from  a  successful  forecast  by  RLAM  during 
February,  1979: 

(Tk-Ck)  =  (-4.2932  X  10‘2,  -  9.8038  X  10'2,  -.105177,  -4.1259  x  10  ", 

2.4924  X  10‘2,  7.1272  X  10'3,  2.1306  X  10  2, 

-7.9304  X  10'2,  -2.7754  X  10'2,  -.  17636,  .  38922) 

If  we  assume,  as  BMRC  did,  that  Ao  =  ^  and  ok  =  ^(dk 
+  crk<1),  and  that  W a'1  is  given  by  averaging  level  W  and 
dividing  by  layer  a  one  obtains 

(Wk)  =  (-6.038  X  10'7,  -9.967  X  10'7,  -2.340  X  10'6,  -8.988  X  10  \ 

-2.381  X  10'7,  -1.389  X  10'7,  -7.609  x  10  7 , 

-1.488  X  10'6,  6.309  x  10  7) 

The  first  value  is  the  tendency  of  In  p  and  is  of 
acceptable  magnitude.  When  the  averaged  d  were  used, 
however,  (Wk)  became 

(-1.518  x  10'5,  1.238  x  10‘5,  -1.444  x  10'5,  9.923  x  10'6, 

9.223  X  10  6,  -7.421  X  10'6,  6.300  X  10‘6,  -4.763  x  10  6, 

2.801  X  10'6,  -3.667  x  10'6,  1.262  X  10'6) 

This  increase  by  two  orders  of  magnitude  of  the  surface 
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pressure  tendency  yields  values  too  large  for  stability  and 

graphically  portrays  the  sensitivity  of  inverting  the  ^ 

matrix.  The  same  occurred  when  the  NMC  12-layer  structure 

was  assumed.  Here  La  was  highly  resolved  near  the  surface, 

less  resolved  in  the  mid-troposhpere  and  highly  resolved  ' 

near  the  tropopause.  The  relationship  of  a  to  a  is  o£  = 

(a/1*0  -  ak+,<U°)/[  (3+k)  (ak  -  dk.,)].  When  a  was 

used  WT  was  -3.896  x  10'7,  while  averaged  a  values  resulted 
in  -1.518  x  10'5,  again  a  two  order-of-magnitude  increase, 
albeit  not  as  large  as  with  the  BMRC  discretization. 

5.  Conclusion 

Apparently,  the  BMRC  scheme  can  lead  to  instabilities 
wrought  by  the  vertical  differencing  and  the  necessity  to 
derive  pressure  tendencies  from  the  vertical  profile  of  the 
temperature  tendency  minus  the  non-linear  terms  of  the 
thermodynamic  equation.  This  linkage,  although 
theoretically  sound,  depends  inordinately  on  the 

differencing  scheme,  demonstrating  that  the  fa  matrix  may  be 

ill-conditioned  and  unreliable  as  a  means  for  deriving 

surface  pressure  and  vertical  velocity.  Inverting  the 

problem  is  a  possibility.  It  involves  substituting  in  (5)  *■ 

for  <p  in  terms  of  W  by  means  of  (7)  and  (9)  and  solving  for 

—  T  — f 

W  then  deriving  T  .  Unfortunately,  that  procedure  does  not 
solve  the  problem,  first,  because  (9)  is  still  used  in  the 
substitution  and  second,  because  values  for  W  at  the  lateral 
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boundaries  are  hard  to  come  by  since  vertical  velocity  is 
normally  derived  diagnostically  from  forecast  divergence  and 
is  not  stored  as  a  prognostic  variable.  Attempts  at 

,  — f  _T 

forecasting  W  and  deriving  T  did  not  seem  any  more 
satisfactory. 

Despite  its  failings,  the  BMRC  semi-implicit  scheme 
does  have  its  application.  If  time  steps  are  constrained, 
the  scheme  could  produce  forecasts  still  more  economically 
than  explicit  schemes  that  must  keep  their  time  steps  in 
proportion  to  the  grid  size.  The  problem  is  knowing  ahead 
of  time  how  large  a  time  step  the  scheme  can  tolerate. 
Several  RLAM  forecasts  were  produced  over  North  America  and 
the  North  Sea  with  the  semi-implicit  scheme.  These 
forecasts  contained  physical  parameterizations  including 
boundary  layer  fluxes  and  large-scale  and  convective 
precipitation  with  a  resolution  of  approximately  200  km  over 
North  America  and  about  100  km  over  the  North  Sea.  The 
semi-implicit  scheme  yielded  forecasts  with  time  steps  of 
about  600  s  for  both  locations  and  these  forecasts  compared 
well  with  the  Brown-Campana  explicit  formulation  using  time 
steps  of  180  s  for  North  America  and  120  s  for  the  North 
Sea.  Near  the  lateral  boundaries,  the  North  American 
forecasts  seemed  to  be  noisier  than  the  Brown-Campana 
forecasts.  But  even  with  a  time  step  of  600  s,  one  forecast 
over  the  North  Sea  failed  after  successfully  forecasting  for 
close  to  40  h  of  simulated  time.  When  the  forecast  was  re- 
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run  with  a  time  step  of  480  s,  it  was  stable  to  48  h.  This 
emphasizes  the  unreliability  of  the  scheme  and  serves  as  an 
illustration  of  what  can  go  wrong  without  warning  because  of 
an  unexpected  shock  to  the  model. 

The  failure  of  the  BMRC  scheme  to  provide  stable 
results  is  not  necessarily  the  result  of  a  flaw  endemic  to 
all  semi-implicit  schemes.  As  mentioned,  the  semi-implicit 
scheme  of  the  AFGL  gsm  or  of  the  Canadian  model  are 
apparently  not  subject  to  the  same  difficulties  as  the  BMRC 
model.  It  may  even  be  possible  to  slightly  modify  the  BMRC 
scheme  to  make  it  less  sensitive.  Results  from  this  study, 
however,  should  encourage  close  scrutiny  of  any  numerical 
scheme  before  adaptation  to  any  model . 
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